Model Expansion
Alex Kale
2021-02-01
Model expansion: Building out a linear in log odds model
Load data.
df <- read_csv("pilot-anonymous.csv")## Parsed with column specification:
## cols(
## .default = col_double(),
## workerId = col_character(),
## condition = col_character(),
## interact = col_logical(),
## gender = col_character(),
## age = col_character(),
## education = col_character(),
## chart_use = col_character(),
## problems = col_character(),
## interactions = col_character(),
## trial = col_character()
## )## See spec(...) for full column specifications.head(df)## # A tibble: 6 x 34
## workerId batch bonus condition duration interact n_data_conds n_trials gender
## <chr> <dbl> <dbl> <chr> <dbl> <lgl> <dbl> <dbl> <chr>
## 1 fccb21d5 1 0.25 icons 1269. FALSE 16 16 man
## 2 fccb21d5 1 0.25 icons 1269. FALSE 16 16 man
## 3 fccb21d5 1 0.25 icons 1269. FALSE 16 16 man
## 4 fccb21d5 1 0.25 icons 1269. FALSE 16 16 man
## 5 fccb21d5 1 0.25 icons 1269. FALSE 16 16 man
## 6 fccb21d5 1 0.25 icons 1269. FALSE 16 16 man
## # … with 25 more variables: age <chr>, education <chr>, chart_use <chr>,
## # problems <chr>, A <dbl>, B <dbl>, C <dbl>, D <dbl>, abs_err <dbl>,
## # causal_support <dbl>, ground_truth <dbl>, interactions <chr>, n <dbl>,
## # nA <dbl>, nB <dbl>, nC <dbl>, nD <dbl>, p_immun <dbl>, payoff <dbl>,
## # q_idx <dbl>, response_A <dbl>, response_B <dbl>, trial <chr>,
## # trial_dur <dbl>, trial_idx <dbl>Calculate a log response ratio lrr to model as a function of causal_support. Also, convert predictor variables to factors for modeling if need be.
model_df <- df %>%
unite("vis", condition, interact, sep = "_") %>%
mutate(
# response units
response_A = if_else(
response_A > 99.5, 99.5,
if_else(
response_A < 0.5, 0.5,
as.numeric(response_A))),
response_B = if_else(
response_B > 99.5, 99.5,
if_else(
response_B < 0.5, 0.5,
as.numeric(response_B))),
lrr = log(response_A / 100) - log(response_B / 100),
# predictors as factors
worker = as.factor(workerId),
vis = as.factor(vis),
n = as.factor(n),
# derived predictors
delta_p = (C + D)/(nC + nD) - (A + B)/(nA + nB)
) Simple slope and intercept
Prior predictive check
p1 <- brm(data = model_df, family = "gaussian",
lrr ~ causal_support,
prior = c(prior(normal(-0.21, 1), class = Intercept), # center at qlogis(mean(model_df$response_A) / 100)
prior(normal(1, 2), class = b, coef = causal_support), # center at unbiased slope
prior(normal(0, 1), class = sigma)), # weakly informative half-normal
sample_prior = "only",
iter = 3000, warmup = 500, chains = 2, cores = 2)## Compiling the C++ model## Start samplingmodel_df %>%
select(causal_support, response_A) %>%
add_predicted_draws(p1, prediction = "lrr_rep", seed = 1234, n = 500) %>%
mutate(
# transform to probability units
response_A_rep = plogis(lrr_rep) * 100
) %>%
ggplot(aes(x = causal_support, y = response_A_rep)) +
stat_lineribbon(.width = c(0.95, 0.8, 0.5)) +
geom_point(data = model_df, aes(y = response_A), alpha = 0.5) +
scale_fill_brewer() +
labs(subtitle = "Prior predictive distribution") +
theme(panel.grid = element_blank())
Fit model
m1 <- brm(data = model_df, family = "gaussian",
lrr ~ causal_support,
prior = c(prior(normal(-0.21, 1), class = Intercept), # center at qlogis(mean(model_df$response_A) / 100)
prior(normal(1, 2), class = b, coef = causal_support), # center at unbiased slope
prior(normal(0, 1), class = sigma)), # weakly informative half-normal
iter = 3000, warmup = 500, chains = 2, cores = 2,
file = "model-fits/1_simple")Diagnostics
summary(m1)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: lrr ~ causal_support
## Data: model_df (Number of observations: 306)
## Samples: 2 chains, each with iter = 3000; warmup = 500; thin = 1;
## total post-warmup samples = 5000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -0.62 0.10 -0.81 -0.43 1.00 5337 3011
## causal_support 0.08 0.01 0.06 0.11 1.00 6992 4198
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.61 0.06 1.49 1.74 1.00 4959 4021
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m1)
pairs(m1)
Posterior predictive
model_df %>%
select(causal_support, response_A) %>%
add_predicted_draws(m1, prediction = "lrr_rep", seed = 1234, n = 500) %>%
mutate(
# transform to probability units
response_A_rep = plogis(lrr_rep) * 100
) %>%
ggplot(aes(x = causal_support, y = response_A_rep)) +
stat_lineribbon(.width = c(0.95, 0.8, 0.5)) +
geom_point(data = model_df, aes(y = response_A), alpha = 0.5) +
scale_fill_brewer() +
labs(subtitle = "Posterior predictive distribution") +
theme(panel.grid = element_blank())
Adding effects of data manipulations
Prior predictive check
p2 <- brm(data = model_df, family = "gaussian",
lrr ~ causal_support*delta_p*n,
prior = c(prior(normal(-0.21, 1), class = Intercept), # center at qlogis(mean(model_df$response_A) / 100)
prior(normal(0, 1), class = b), # center predictor effects at 0
prior(normal(1, 2), class = b, coef = causal_support), # center at unbiased slope
prior(normal(0, 1), class = sigma)), # weakly informative half-normal
sample_prior = "only",
iter = 3000, warmup = 500, chains = 2, cores = 2)## Compiling the C++ model## Start samplingexpand_grid(
causal_support = quantile(model_df$causal_support, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 50))),
delta_p = quantile(model_df$delta_p, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 20))),
n = unique(model_df$n)) %>%
add_predicted_draws(p2, prediction = "lrr_rep", seed = 1234, n = 500) %>%
mutate(
# transform to probability units
response_A_rep = plogis(lrr_rep) * 100
) %>%
ggplot(aes(x = causal_support, y = response_A_rep)) +
stat_lineribbon(.width = c(0.95, 0.8, 0.5)) +
geom_point(data = model_df, aes(y = response_A), alpha = 0.5) +
scale_fill_brewer() +
labs(subtitle = "Prior predictive distribution") +
theme(panel.grid = element_blank()) +
facet_grid(n ~ .)
Fit model
m2 <- brm(data = model_df, family = "gaussian",
lrr ~ causal_support*delta_p*n,
prior = c(prior(normal(-0.21, 1), class = Intercept), # center at qlogis(mean(model_df$response_A) / 100)
prior(normal(0, 1), class = b), # center predictor effects at 0
prior(normal(1, 2), class = b, coef = causal_support), # center at unbiased slope
prior(normal(0, 1), class = sigma)), # weakly informative half-normal
iter = 3000, warmup = 500, chains = 2, cores = 2,
file = "model-fits/2_data-conds")Diagnostics
summary(m2)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: lrr ~ causal_support * delta_p * n
## Data: model_df (Number of observations: 306)
## Samples: 2 chains, each with iter = 3000; warmup = 500; thin = 1;
## total post-warmup samples = 5000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -0.78 0.18 -1.13 -0.43 1.00 3655
## causal_support 0.40 0.17 0.06 0.72 1.00 1396
## delta_p 0.36 0.97 -1.57 2.19 1.00 5381
## n500 0.28 0.25 -0.22 0.76 1.00 3730
## n1000 0.21 0.26 -0.31 0.70 1.00 3771
## n5000 0.30 0.26 -0.22 0.82 1.00 4468
## causal_support:delta_p -0.93 0.68 -2.24 0.42 1.00 4779
## causal_support:n500 -0.12 0.19 -0.50 0.26 1.00 1633
## causal_support:n1000 0.00 0.18 -0.34 0.36 1.00 1511
## causal_support:n5000 -0.31 0.17 -0.64 0.03 1.00 1386
## delta_p:n500 0.01 0.97 -1.87 1.90 1.00 5595
## delta_p:n1000 0.06 0.99 -1.87 2.02 1.00 8708
## delta_p:n5000 0.10 1.02 -1.92 2.08 1.00 6118
## causal_support:delta_p:n500 -0.11 0.96 -1.95 1.77 1.00 5464
## causal_support:delta_p:n1000 -1.01 0.91 -2.80 0.77 1.00 5441
## causal_support:delta_p:n5000 0.39 0.75 -1.06 1.85 1.00 4459
## Tail_ESS
## Intercept 3613
## causal_support 2283
## delta_p 3262
## n500 3434
## n1000 3714
## n5000 3908
## causal_support:delta_p 3745
## causal_support:n500 2651
## causal_support:n1000 2493
## causal_support:n5000 2170
## delta_p:n500 3811
## delta_p:n1000 4044
## delta_p:n5000 3733
## causal_support:delta_p:n500 3211
## causal_support:delta_p:n1000 3985
## causal_support:delta_p:n5000 4070
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.55 0.06 1.44 1.68 1.00 6169 3856
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m2)
# intercepts
pairs(m2, pars = c("b_Intercept", "b_delta_p", "b_n500", "b_n1000", "b_n5000"))
# slopes on causal support
pairs(m2, pars = c("b_causal_support", "b_causal_support:delta_p", "b_causal_support:n500", "b_causal_support:n1000", "b_causal_support:n5000"))
# slopes on delta_p and sigma
pairs(m2, pars = c("b_delta_p", "b_delta_p:n500", "b_delta_p:n1000", "b_delta_p:n5000", "sigma"))
# slope interactions
pairs(m2, pars = c("b_causal_support:delta_p", "b_causal_support:delta_p:n500", "b_causal_support:delta_p:n1000", "b_causal_support:delta_p:n5000"))
Posterior predictive
expand_grid(
causal_support = quantile(model_df$causal_support, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 50))),
delta_p = quantile(model_df$delta_p, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 20))),
n = unique(model_df$n)) %>%
add_predicted_draws(m2, prediction = "lrr_rep", seed = 1234, n = 500) %>%
mutate(
# transform to probability units
response_A_rep = plogis(lrr_rep) * 100
) %>%
ggplot(aes(x = causal_support, y = response_A_rep)) +
stat_lineribbon(.width = c(0.95, 0.8, 0.5)) +
geom_point(data = model_df, aes(y = response_A), alpha = 0.5) +
scale_fill_brewer() +
labs(subtitle = "Posterior predictive distribution") +
theme(panel.grid = element_blank()) +
facet_grid(n ~ .)
Adding effects of visualization conditions
Prior predictive check
p3 <- brm(data = model_df, family = "gaussian",
lrr ~ causal_support*delta_p*n*vis,
prior = c(prior(normal(-0.21, 1), class = Intercept), # center at qlogis(mean(model_df$response_A) / 100)
prior(normal(0, 1), class = b), # center predictor effects at 0
prior(normal(1, 2), class = b, coef = causal_support), # center at unbiased slope
prior(normal(0, 1), class = sigma)), # weakly informative half-normal
sample_prior = "only",
iter = 3000, warmup = 500, chains = 2, cores = 2)## Compiling the C++ model## Start samplingexpand_grid(
causal_support = quantile(model_df$causal_support, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 50))),
delta_p = quantile(model_df$delta_p, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 20))),
n = unique(model_df$n),
vis = unique(model_df$vis)) %>%
add_predicted_draws(p3, prediction = "lrr_rep", seed = 1234, n = 500) %>%
mutate(
# transform to probability units
response_A_rep = plogis(lrr_rep) * 100
) %>%
ggplot(aes(x = causal_support, y = response_A_rep)) +
stat_lineribbon(.width = c(0.95, 0.8, 0.5)) +
geom_point(data = model_df, aes(y = response_A), alpha = 0.5) +
scale_fill_brewer() +
labs(subtitle = "Prior predictive distribution") +
theme(panel.grid = element_blank()) +
facet_grid(n ~ vis)
Fit model
m3 <- brm(data = model_df, family = "gaussian",
lrr ~ causal_support*delta_p*n*vis,
prior = c(prior(normal(-0.21, 1), class = Intercept), # center at qlogis(mean(model_df$response_A) / 100)
prior(normal(0, 1), class = b), # center predictor effects at 0
prior(normal(1, 2), class = b, coef = causal_support), # center at unbiased slope
prior(normal(0, 1), class = sigma)), # weakly informative half-normal
iter = 3000, warmup = 500, chains = 2, cores = 2,
file = "model-fits/3_vis")Diagnostics
summary(m3)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: lrr ~ causal_support * delta_p * n * vis
## Data: model_df (Number of observations: 306)
## Samples: 2 chains, each with iter = 3000; warmup = 500; thin = 1;
## total post-warmup samples = 5000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## Intercept -0.57 0.22 -1.00 -0.14
## causal_support 0.51 0.20 0.13 0.91
## delta_p 0.41 1.00 -1.57 2.39
## n500 0.07 0.31 -0.55 0.66
## n1000 -0.21 0.32 -0.85 0.40
## n5000 0.42 0.32 -0.19 1.04
## vistext_FALSE -0.46 0.29 -1.03 0.11
## causal_support:delta_p -1.03 0.73 -2.48 0.45
## causal_support:n500 -0.21 0.23 -0.67 0.24
## causal_support:n1000 -0.09 0.21 -0.52 0.32
## causal_support:n5000 -0.44 0.20 -0.84 -0.06
## delta_p:n500 0.01 0.99 -1.94 1.88
## delta_p:n1000 0.07 1.03 -1.96 2.02
## delta_p:n5000 0.10 1.00 -1.87 2.03
## causal_support:vistext_FALSE -0.28 0.27 -0.81 0.25
## delta_p:vistext_FALSE -0.09 0.98 -2.06 1.78
## n500:vistext_FALSE 0.48 0.41 -0.30 1.29
## n1000:vistext_FALSE 0.94 0.43 0.10 1.78
## n5000:vistext_FALSE -0.22 0.43 -1.07 0.62
## causal_support:delta_p:n500 -0.12 0.94 -1.96 1.70
## causal_support:delta_p:n1000 -1.01 0.90 -2.75 0.82
## causal_support:delta_p:n5000 0.20 0.78 -1.32 1.72
## causal_support:delta_p:vistext_FALSE -0.12 0.82 -1.74 1.47
## causal_support:n500:vistext_FALSE 0.25 0.31 -0.36 0.87
## causal_support:n1000:vistext_FALSE 0.28 0.29 -0.27 0.85
## causal_support:n5000:vistext_FALSE 0.35 0.27 -0.18 0.88
## delta_p:n500:vistext_FALSE -0.08 1.01 -2.06 1.84
## delta_p:n1000:vistext_FALSE -0.02 0.98 -1.91 1.90
## delta_p:n5000:vistext_FALSE 0.01 1.02 -1.98 2.00
## causal_support:delta_p:n500:vistext_FALSE 0.23 0.97 -1.64 2.14
## causal_support:delta_p:n1000:vistext_FALSE -0.23 0.97 -2.20 1.64
## causal_support:delta_p:n5000:vistext_FALSE -0.04 0.86 -1.72 1.69
## Rhat Bulk_ESS Tail_ESS
## Intercept 1.00 3519 4156
## causal_support 1.00 1402 2427
## delta_p 1.00 8551 3954
## n500 1.00 3897 3925
## n1000 1.00 3865 3775
## n5000 1.00 4271 4117
## vistext_FALSE 1.00 3103 3584
## causal_support:delta_p 1.00 5276 4038
## causal_support:n500 1.00 1735 2940
## causal_support:n1000 1.00 1559 2623
## causal_support:n5000 1.00 1412 2479
## delta_p:n500 1.00 6487 3547
## delta_p:n1000 1.00 7262 3631
## delta_p:n5000 1.00 6861 3882
## causal_support:vistext_FALSE 1.00 1429 2666
## delta_p:vistext_FALSE 1.00 7372 3700
## n500:vistext_FALSE 1.00 3458 4038
## n1000:vistext_FALSE 1.00 3913 3647
## n5000:vistext_FALSE 1.00 3957 3950
## causal_support:delta_p:n500 1.00 7322 3461
## causal_support:delta_p:n1000 1.00 6163 4096
## causal_support:delta_p:n5000 1.00 5192 4015
## causal_support:delta_p:vistext_FALSE 1.00 8452 4146
## causal_support:n500:vistext_FALSE 1.00 1814 2937
## causal_support:n1000:vistext_FALSE 1.00 1539 2759
## causal_support:n5000:vistext_FALSE 1.00 1447 2552
## delta_p:n500:vistext_FALSE 1.00 7424 3996
## delta_p:n1000:vistext_FALSE 1.00 7829 3667
## delta_p:n5000:vistext_FALSE 1.00 7774 3986
## causal_support:delta_p:n500:vistext_FALSE 1.00 7063 3514
## causal_support:delta_p:n1000:vistext_FALSE 1.00 8150 3422
## causal_support:delta_p:n5000:vistext_FALSE 1.00 6815 4045
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.54 0.07 1.41 1.67 1.00 7060 3786
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m3)
Too many pairs to plot, so we’ll only do it as needed to check issues from here on out.
Posterior predictive
expand_grid(
causal_support = quantile(model_df$causal_support, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 20))),
delta_p = quantile(model_df$delta_p, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 20))),
n = unique(model_df$n),
vis = unique(model_df$vis)) %>%
add_predicted_draws(m3, prediction = "lrr_rep", seed = 1234, n = 500) %>%
mutate(
# transform to probability units
response_A_rep = plogis(lrr_rep) * 100
) %>%
ggplot(aes(x = causal_support, y = response_A_rep)) +
stat_lineribbon(.width = c(0.95, 0.8, 0.5)) +
geom_point(data = model_df, aes(y = response_A), alpha = 0.5) +
scale_fill_brewer() +
labs(subtitle = "Posterior predictive distribution") +
theme(panel.grid = element_blank()) +
facet_grid(n ~ vis)
Adding variance effects
Prior predictive check
p4 <- brm(data = model_df, family = "gaussian",
formula = bf(lrr ~ causal_support*delta_p*n*vis,
sigma ~ abs(causal_support)),
prior = c(prior(normal(-0.21, 1), class = Intercept), # center at qlogis(mean(model_df$response_A) / 100)
prior(normal(0, 1), class = b), # center predictor effects at 0
prior(normal(1, 2), class = b, coef = causal_support), # center at unbiased slope
prior(normal(0, 1), class = b, dpar = sigma), # weakly informative half-normal
prior(normal(0, 1), class = Intercept, dpar = sigma)), # weakly informative half-normal
sample_prior = "only",
iter = 3000, warmup = 500, chains = 2, cores = 2)## Compiling the C++ model## Start samplingexpand_grid(
causal_support = quantile(model_df$causal_support, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 50))),
delta_p = quantile(model_df$delta_p, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 20))),
n = unique(model_df$n),
vis = unique(model_df$vis)) %>%
add_predicted_draws(p4, prediction = "lrr_rep", seed = 1234, n = 500) %>%
mutate(
# transform to probability units
response_A_rep = plogis(lrr_rep) * 100
) %>%
ggplot(aes(x = causal_support, y = response_A_rep)) +
stat_lineribbon(.width = c(0.95, 0.8, 0.5)) +
geom_point(data = model_df, aes(y = response_A), alpha = 0.5) +
scale_fill_brewer() +
labs(subtitle = "Prior predictive distribution") +
theme(panel.grid = element_blank()) +
facet_grid(n ~ vis)
Fit model
m4 <- brm(data = model_df, family = "gaussian",
formula = bf(lrr ~ causal_support*delta_p*n*vis,
sigma ~ abs(causal_support)),
prior = c(prior(normal(-0.21, 1), class = Intercept), # center at qlogis(mean(model_df$response_A) / 100)
prior(normal(0, 1), class = b), # center predictor effects at 0
prior(normal(1, 2), class = b, coef = causal_support), # center at unbiased slope
prior(normal(0, 1), class = b, dpar = sigma), # weakly informative half-normal
prior(normal(0, 1), class = Intercept, dpar = sigma)), # weakly informative half-normal
iter = 3000, warmup = 500, chains = 2, cores = 2,
file = "model-fits/4_var")Diagnostics
summary(m4)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: lrr ~ causal_support * delta_p * n * vis
## sigma ~ abs(causal_support)
## Data: model_df (Number of observations: 306)
## Samples: 2 chains, each with iter = 3000; warmup = 500; thin = 1;
## total post-warmup samples = 5000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## Intercept -0.58 0.23 -1.03 -0.13
## sigma_Intercept 0.46 0.05 0.37 0.56
## causal_support 0.52 0.20 0.12 0.92
## delta_p 0.37 0.95 -1.50 2.24
## n500 0.07 0.32 -0.54 0.70
## n1000 -0.21 0.33 -0.86 0.45
## n5000 0.47 0.32 -0.16 1.12
## vistext_FALSE -0.45 0.31 -1.06 0.16
## causal_support:delta_p -1.06 0.72 -2.45 0.36
## causal_support:n500 -0.22 0.24 -0.70 0.25
## causal_support:n1000 -0.10 0.22 -0.51 0.33
## causal_support:n5000 -0.46 0.20 -0.86 -0.05
## delta_p:n500 -0.00 1.01 -2.00 1.97
## delta_p:n1000 0.04 1.01 -1.91 2.00
## delta_p:n5000 0.09 0.98 -1.84 2.01
## causal_support:vistext_FALSE -0.29 0.28 -0.82 0.25
## delta_p:vistext_FALSE -0.09 0.99 -2.09 1.82
## n500:vistext_FALSE 0.48 0.43 -0.36 1.33
## n1000:vistext_FALSE 0.91 0.44 0.05 1.79
## n5000:vistext_FALSE -0.32 0.45 -1.21 0.57
## causal_support:delta_p:n500 -0.13 0.96 -2.00 1.76
## causal_support:delta_p:n1000 -0.99 0.91 -2.73 0.78
## causal_support:delta_p:n5000 0.17 0.75 -1.31 1.66
## causal_support:delta_p:vistext_FALSE -0.14 0.84 -1.81 1.49
## causal_support:n500:vistext_FALSE 0.26 0.32 -0.36 0.90
## causal_support:n1000:vistext_FALSE 0.29 0.29 -0.27 0.86
## causal_support:n5000:vistext_FALSE 0.36 0.28 -0.17 0.89
## delta_p:n500:vistext_FALSE -0.04 0.99 -2.02 1.91
## delta_p:n1000:vistext_FALSE -0.01 1.02 -2.03 2.04
## delta_p:n5000:vistext_FALSE 0.04 1.02 -1.95 2.01
## causal_support:delta_p:n500:vistext_FALSE 0.23 0.99 -1.74 2.17
## causal_support:delta_p:n1000:vistext_FALSE -0.20 0.99 -2.10 1.75
## causal_support:delta_p:n5000:vistext_FALSE -0.04 0.86 -1.72 1.65
## sigma_abscausal_support -0.01 0.01 -0.02 0.00
## Rhat Bulk_ESS Tail_ESS
## Intercept 1.00 2882 3040
## sigma_Intercept 1.00 5472 3979
## causal_support 1.00 1215 1850
## delta_p 1.00 5249 3460
## n500 1.00 3159 3451
## n1000 1.00 2672 3183
## n5000 1.00 3049 3326
## vistext_FALSE 1.00 2433 3057
## causal_support:delta_p 1.00 4032 3966
## causal_support:n500 1.00 1526 2239
## causal_support:n1000 1.00 1303 2108
## causal_support:n5000 1.00 1215 1803
## delta_p:n500 1.00 6844 3544
## delta_p:n1000 1.00 5725 3444
## delta_p:n5000 1.00 5867 3867
## causal_support:vistext_FALSE 1.00 1415 2465
## delta_p:vistext_FALSE 1.00 4855 3689
## n500:vistext_FALSE 1.00 2857 3367
## n1000:vistext_FALSE 1.00 2543 3497
## n5000:vistext_FALSE 1.00 3013 3366
## causal_support:delta_p:n500 1.00 7162 3377
## causal_support:delta_p:n1000 1.00 4008 3544
## causal_support:delta_p:n5000 1.00 4383 3936
## causal_support:delta_p:vistext_FALSE 1.00 4328 3516
## causal_support:n500:vistext_FALSE 1.00 1693 2761
## causal_support:n1000:vistext_FALSE 1.00 1538 2503
## causal_support:n5000:vistext_FALSE 1.00 1439 2440
## delta_p:n500:vistext_FALSE 1.00 5331 3691
## delta_p:n1000:vistext_FALSE 1.00 6560 3663
## delta_p:n5000:vistext_FALSE 1.00 7370 3901
## causal_support:delta_p:n500:vistext_FALSE 1.00 5940 3678
## causal_support:delta_p:n1000:vistext_FALSE 1.00 5921 3474
## causal_support:delta_p:n5000:vistext_FALSE 1.00 5550 3216
## sigma_abscausal_support 1.00 8307 3729
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m4)
Posterior predictive
expand_grid(
causal_support = quantile(model_df$causal_support, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 20))),
delta_p = quantile(model_df$delta_p, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 20))),
n = unique(model_df$n),
vis = unique(model_df$vis)) %>%
add_predicted_draws(m4, prediction = "lrr_rep", seed = 1234, n = 500) %>%
mutate(
# transform to probability units
response_A_rep = plogis(lrr_rep) * 100
) %>%
ggplot(aes(x = causal_support, y = response_A_rep)) +
stat_lineribbon(.width = c(0.95, 0.8, 0.5)) +
geom_point(data = model_df, aes(y = response_A), alpha = 0.5) +
scale_fill_brewer() +
labs(subtitle = "Posterior predictive distribution") +
theme(panel.grid = element_blank()) +
facet_grid(n ~ vis)
Adding random effects on slope and intercept
Prior predictive check. Now that we’ve added hierarchy we’ll narrow our priors to get a bit of regularization.
p5 <- brm(data = model_df, family = "gaussian",
formula = bf(lrr ~ causal_support*delta_p*n*vis + (causal_support|workerId),
sigma ~ abs(causal_support) + (1|workerId)),# + (abs(causal_support)|workerId)),
prior = c(prior(normal(-0.21, 1), class = Intercept), # center at qlogis(mean(model_df$response_A) / 100)
prior(normal(0, 0.5), class = b), # center predictor effects at 0
prior(normal(1, 0.5), class = b, coef = causal_support), # center at unbiased slope
prior(normal(0, 0.2), class = b, dpar = sigma), # weakly informative half-normal
prior(normal(0, 0.2), class = Intercept, dpar = sigma), # weakly informative half-normal
prior(normal(0, 0.2), class = sd), # weakly informative half-normal
prior(lkj(4), class = cor)), # avoiding large correlations
sample_prior = "only",
iter = 3000, warmup = 500, chains = 2, cores = 2)## Compiling the C++ model## Start sampling## Warning: There were 2 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup## Warning: Examine the pairs() plot to diagnose sampling problemsmodel_df %>%
group_by(n, vis, workerId) %>%
data_grid(
causal_support = quantile(model_df$causal_support, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 20))),
delta_p = quantile(model_df$delta_p, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 10)))) %>%
add_predicted_draws(p5, prediction = "lrr_rep", seed = 1234, n = 500) %>%
mutate(
# transform to probability units
response_A_rep = plogis(lrr_rep) * 100
) %>%
ggplot(aes(x = causal_support, y = response_A_rep, fill = vis)) +
stat_lineribbon(.width = c(0.95)) +
geom_point(data = model_df, aes(y = response_A), alpha = 0.5) +
scale_fill_brewer(type = "qual", palette = 2) +
scale_color_brewer(type = "qual", palette = 2) +
labs(subtitle = "Prior predictive distribution") +
theme(panel.grid = element_blank()) +
facet_wrap(workerId ~ .)
Fit model. Random slope effects on sigma did not work out.
m5 <- brm(data = model_df, family = "gaussian",
formula = bf(lrr ~ causal_support*delta_p*n*vis + (causal_support|workerId),
sigma ~ abs(causal_support) + (1|workerId)),# + (abs(causal_support)|workerId)),
prior = c(prior(normal(-0.21, 1), class = Intercept), # center at qlogis(mean(model_df$response_A) / 100)
prior(normal(0, 0.5), class = b), # center predictor effects at 0
prior(normal(1, 0.5), class = b, coef = causal_support), # center at unbiased slope
prior(normal(0, 0.2), class = b, dpar = sigma), # weakly informative half-normal
prior(normal(0, 0.2), class = Intercept, dpar = sigma), # weakly informative half-normal
prior(normal(0, 0.2), class = sd), # weakly informative half-normal
prior(lkj(4), class = cor)), # avoiding large correlations
iter = 3000, warmup = 500, chains = 2, cores = 2,
control = list(adapt_delta = 0.99, max_treedepth = 12),
file = "model-fits/5b_re-simple")Diagnostics
summary(m5)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: lrr ~ causal_support * delta_p * n * vis + (causal_support | workerId)
## sigma ~ abs(causal_support) + (1 | workerId)
## Data: model_df (Number of observations: 306)
## Samples: 2 chains, each with iter = 3000; warmup = 500; thin = 1;
## total post-warmup samples = 5000
##
## Group-Level Effects:
## ~workerId (Number of levels: 18)
## Estimate Est.Error l-95% CI u-95% CI Rhat
## sd(Intercept) 0.63 0.09 0.46 0.80 1.00
## sd(causal_support) 0.07 0.02 0.04 0.12 1.00
## sd(sigma_Intercept) 0.46 0.10 0.29 0.70 1.00
## cor(Intercept,causal_support) -0.43 0.17 -0.73 -0.05 1.00
## Bulk_ESS Tail_ESS
## sd(Intercept) 2836 3542
## sd(causal_support) 1896 3050
## sd(sigma_Intercept) 1516 2534
## cor(Intercept,causal_support) 2518 3213
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## Intercept -0.58 0.24 -1.07 -0.11
## sigma_Intercept 0.17 0.11 -0.05 0.38
## causal_support 0.34 0.14 0.07 0.63
## delta_p 0.17 0.50 -0.79 1.15
## n500 -0.04 0.18 -0.40 0.31
## n1000 0.10 0.18 -0.28 0.45
## n5000 0.28 0.19 -0.09 0.65
## vistext_FALSE -0.19 0.31 -0.77 0.42
## causal_support:delta_p -0.60 0.39 -1.35 0.17
## causal_support:n500 -0.18 0.15 -0.48 0.12
## causal_support:n1000 -0.14 0.15 -0.44 0.14
## causal_support:n5000 -0.27 0.14 -0.55 0.00
## delta_p:n500 0.09 0.50 -0.87 1.06
## delta_p:n1000 0.01 0.50 -0.94 0.94
## delta_p:n5000 0.04 0.51 -0.94 1.02
## causal_support:vistext_FALSE 0.01 0.21 -0.41 0.40
## delta_p:vistext_FALSE 0.01 0.51 -0.98 1.05
## n500:vistext_FALSE 0.20 0.29 -0.37 0.76
## n1000:vistext_FALSE 0.27 0.29 -0.29 0.84
## n5000:vistext_FALSE -0.06 0.29 -0.63 0.52
## causal_support:delta_p:n500 -0.05 0.47 -0.93 0.89
## causal_support:delta_p:n1000 -0.40 0.47 -1.30 0.51
## causal_support:delta_p:n5000 0.06 0.42 -0.77 0.85
## causal_support:delta_p:vistext_FALSE -0.20 0.46 -1.09 0.71
## causal_support:n500:vistext_FALSE 0.18 0.22 -0.26 0.63
## causal_support:n1000:vistext_FALSE 0.19 0.21 -0.22 0.60
## causal_support:n5000:vistext_FALSE 0.03 0.20 -0.37 0.43
## delta_p:n500:vistext_FALSE 0.01 0.50 -0.97 1.01
## delta_p:n1000:vistext_FALSE -0.00 0.49 -0.97 0.96
## delta_p:n5000:vistext_FALSE 0.01 0.51 -0.98 1.01
## causal_support:delta_p:n500:vistext_FALSE 0.06 0.50 -0.93 1.04
## causal_support:delta_p:n1000:vistext_FALSE -0.10 0.49 -1.07 0.88
## causal_support:delta_p:n5000:vistext_FALSE -0.06 0.47 -0.99 0.84
## sigma_abscausal_support -0.02 0.01 -0.04 0.00
## Rhat Bulk_ESS Tail_ESS
## Intercept 1.00 1771 2867
## sigma_Intercept 1.00 1516 2455
## causal_support 1.00 1479 2454
## delta_p 1.00 6951 3567
## n500 1.00 4332 3753
## n1000 1.00 4318 4220
## n5000 1.00 4300 3605
## vistext_FALSE 1.00 2240 3185
## causal_support:delta_p 1.00 5376 3934
## causal_support:n500 1.00 1671 2895
## causal_support:n1000 1.00 1605 2601
## causal_support:n5000 1.00 1513 2610
## delta_p:n500 1.00 6733 3471
## delta_p:n1000 1.00 7645 3937
## delta_p:n5000 1.00 8993 3933
## causal_support:vistext_FALSE 1.00 1633 2351
## delta_p:vistext_FALSE 1.00 8695 3346
## n500:vistext_FALSE 1.00 4482 4013
## n1000:vistext_FALSE 1.00 4447 3853
## n5000:vistext_FALSE 1.00 3904 3656
## causal_support:delta_p:n500 1.00 6415 3719
## causal_support:delta_p:n1000 1.00 8123 3104
## causal_support:delta_p:n5000 1.00 5846 4157
## causal_support:delta_p:vistext_FALSE 1.00 6713 3904
## causal_support:n500:vistext_FALSE 1.00 1875 2680
## causal_support:n1000:vistext_FALSE 1.00 1799 2443
## causal_support:n5000:vistext_FALSE 1.00 1672 2152
## delta_p:n500:vistext_FALSE 1.00 9810 3611
## delta_p:n1000:vistext_FALSE 1.00 7501 3428
## delta_p:n5000:vistext_FALSE 1.00 9422 3182
## causal_support:delta_p:n500:vistext_FALSE 1.00 9237 3478
## causal_support:delta_p:n1000:vistext_FALSE 1.00 9300 3438
## causal_support:delta_p:n5000:vistext_FALSE 1.00 5302 3693
## sigma_abscausal_support 1.00 3609 3557
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m5)
Posterior predictive
model_df %>%
group_by(n, vis, workerId) %>%
data_grid(
causal_support = quantile(model_df$causal_support, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 20))),
delta_p = quantile(model_df$delta_p, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 10)))) %>%
add_predicted_draws(m5, prediction = "lrr_rep", seed = 1234, n = 500) %>%
mutate(
# transform to probability units
response_A_rep = plogis(lrr_rep) * 100
) %>%
ggplot(aes(x = causal_support, y = response_A_rep, fill = vis)) +
stat_lineribbon(.width = c(0.95)) +
geom_point(data = model_df, aes(y = response_A), alpha = 0.5) +
scale_fill_brewer(type = "qual", palette = 2) +
scale_color_brewer(type = "qual", palette = 2) +
labs(subtitle = "Posterior predictive distribution") +
theme(panel.grid = element_blank()) +
facet_wrap(workerId ~ .)
Adding random effects for interactions of within-subjects manipulations
Same priors as before.
Fit model. Full interaction of random slopes did not work.
m6 <- brm(data = model_df, family = "gaussian",
formula = bf(lrr ~ causal_support*delta_p*n*vis + (causal_support:delta_p + causal_support:n|workerId), #+ (causal_support*delta_p*n|workerId),
sigma ~ abs(causal_support) + (1|workerId)),
prior = c(prior(normal(-0.21, 1), class = Intercept), # center at qlogis(mean(model_df$response_A) / 100)
prior(normal(0, 0.5), class = b), # center predictor effects at 0
prior(normal(1, 0.5), class = b, coef = causal_support), # center at unbiased slope
prior(normal(0, 0.2), class = b, dpar = sigma), # weakly informative half-normal
prior(normal(0, 0.2), class = Intercept, dpar = sigma), # weakly informative half-normal
prior(normal(0, 0.2), class = sd), # weakly informative half-normal
prior(lkj(4), class = cor)), # avoiding large correlations
iter = 3000, warmup = 500, chains = 2, cores = 2,
control = list(adapt_delta = 0.99, max_treedepth = 12),
file = "model-fits/6_re-data_conds")Diagnostics
summary(m6)## Family: gaussian
## Links: mu = identity; sigma = log
## Formula: lrr ~ causal_support * delta_p * n * vis + (causal_support:delta_p + causal_support:n | workerId)
## sigma ~ abs(causal_support) + (1 | workerId)
## Data: model_df (Number of observations: 306)
## Samples: 2 chains, each with iter = 3000; warmup = 500; thin = 1;
## total post-warmup samples = 5000
##
## Group-Level Effects:
## ~workerId (Number of levels: 18)
## Estimate Est.Error l-95% CI
## sd(Intercept) 0.60 0.09 0.45
## sd(causal_support:delta_p) 0.20 0.15 0.01
## sd(causal_support:n500) 0.14 0.09 0.01
## sd(causal_support:n1000) 0.21 0.07 0.08
## sd(causal_support:n5000) 0.05 0.02 0.02
## sd(sigma_Intercept) 0.46 0.11 0.29
## cor(Intercept,causal_support:delta_p) -0.03 0.28 -0.56
## cor(Intercept,causal_support:n500) -0.10 0.25 -0.57
## cor(causal_support:delta_p,causal_support:n500) 0.02 0.29 -0.54
## cor(Intercept,causal_support:n1000) -0.30 0.20 -0.66
## cor(causal_support:delta_p,causal_support:n1000) 0.06 0.29 -0.51
## cor(causal_support:n500,causal_support:n1000) 0.19 0.28 -0.38
## cor(Intercept,causal_support:n5000) -0.35 0.19 -0.69
## cor(causal_support:delta_p,causal_support:n5000) 0.02 0.28 -0.53
## cor(causal_support:n500,causal_support:n5000) 0.17 0.27 -0.41
## cor(causal_support:n1000,causal_support:n5000) 0.41 0.23 -0.09
## u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.79 1.00 3127
## sd(causal_support:delta_p) 0.55 1.00 2070
## sd(causal_support:n500) 0.35 1.00 1328
## sd(causal_support:n1000) 0.36 1.00 1449
## sd(causal_support:n5000) 0.09 1.00 1451
## sd(sigma_Intercept) 0.73 1.00 1526
## cor(Intercept,causal_support:delta_p) 0.53 1.00 5190
## cor(Intercept,causal_support:n500) 0.41 1.00 4120
## cor(causal_support:delta_p,causal_support:n500) 0.56 1.00 3374
## cor(Intercept,causal_support:n1000) 0.12 1.00 3351
## cor(causal_support:delta_p,causal_support:n1000) 0.60 1.00 1857
## cor(causal_support:n500,causal_support:n1000) 0.67 1.00 2152
## cor(Intercept,causal_support:n5000) 0.06 1.00 3074
## cor(causal_support:delta_p,causal_support:n5000) 0.55 1.00 2116
## cor(causal_support:n500,causal_support:n5000) 0.65 1.00 2237
## cor(causal_support:n1000,causal_support:n5000) 0.78 1.00 2971
## Tail_ESS
## sd(Intercept) 2846
## sd(causal_support:delta_p) 2277
## sd(causal_support:n500) 1866
## sd(causal_support:n1000) 1034
## sd(causal_support:n5000) 875
## sd(sigma_Intercept) 2825
## cor(Intercept,causal_support:delta_p) 3431
## cor(Intercept,causal_support:n500) 3457
## cor(causal_support:delta_p,causal_support:n500) 3607
## cor(Intercept,causal_support:n1000) 3368
## cor(causal_support:delta_p,causal_support:n1000) 3036
## cor(causal_support:n500,causal_support:n1000) 3136
## cor(Intercept,causal_support:n5000) 3097
## cor(causal_support:delta_p,causal_support:n5000) 3266
## cor(causal_support:n500,causal_support:n5000) 2970
## cor(causal_support:n1000,causal_support:n5000) 3411
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI
## Intercept -0.56 0.23 -1.01 -0.10
## sigma_Intercept 0.16 0.11 -0.07 0.37
## causal_support 0.35 0.14 0.08 0.62
## delta_p 0.15 0.51 -0.82 1.11
## n500 -0.05 0.18 -0.40 0.31
## n1000 0.12 0.18 -0.24 0.47
## n5000 0.27 0.18 -0.09 0.63
## vistext_FALSE -0.22 0.30 -0.80 0.36
## causal_support:delta_p -0.53 0.41 -1.34 0.27
## causal_support:n500 -0.18 0.16 -0.49 0.15
## causal_support:n1000 -0.06 0.16 -0.37 0.25
## causal_support:n5000 -0.28 0.14 -0.55 -0.01
## delta_p:n500 0.08 0.49 -0.88 1.04
## delta_p:n1000 0.01 0.49 -0.92 0.96
## delta_p:n5000 0.04 0.50 -0.93 1.01
## causal_support:vistext_FALSE -0.00 0.20 -0.39 0.40
## delta_p:vistext_FALSE 0.00 0.50 -0.94 0.99
## n500:vistext_FALSE 0.22 0.29 -0.34 0.77
## n1000:vistext_FALSE 0.29 0.28 -0.26 0.83
## n5000:vistext_FALSE -0.08 0.28 -0.66 0.47
## causal_support:delta_p:n500 -0.07 0.49 -1.00 0.91
## causal_support:delta_p:n1000 -0.25 0.49 -1.21 0.74
## causal_support:delta_p:n5000 0.04 0.41 -0.79 0.85
## causal_support:delta_p:vistext_FALSE -0.14 0.46 -1.03 0.78
## causal_support:n500:vistext_FALSE 0.16 0.23 -0.32 0.61
## causal_support:n1000:vistext_FALSE 0.10 0.23 -0.34 0.54
## causal_support:n5000:vistext_FALSE 0.05 0.21 -0.37 0.44
## delta_p:n500:vistext_FALSE 0.01 0.50 -0.96 0.98
## delta_p:n1000:vistext_FALSE -0.00 0.50 -0.99 0.99
## delta_p:n5000:vistext_FALSE 0.02 0.50 -0.95 1.00
## causal_support:delta_p:n500:vistext_FALSE 0.06 0.50 -0.93 1.04
## causal_support:delta_p:n1000:vistext_FALSE -0.06 0.50 -1.03 0.94
## causal_support:delta_p:n5000:vistext_FALSE -0.06 0.47 -0.97 0.84
## sigma_abscausal_support -0.02 0.01 -0.03 0.00
## Rhat Bulk_ESS Tail_ESS
## Intercept 1.00 1860 2484
## sigma_Intercept 1.00 1591 2126
## causal_support 1.00 1813 2806
## delta_p 1.00 7278 3189
## n500 1.00 4894 3974
## n1000 1.00 4335 4269
## n5000 1.00 3952 3368
## vistext_FALSE 1.00 2550 3360
## causal_support:delta_p 1.00 4834 3546
## causal_support:n500 1.00 2155 3159
## causal_support:n1000 1.00 2056 3182
## causal_support:n5000 1.00 1784 2631
## delta_p:n500 1.00 8715 3410
## delta_p:n1000 1.00 6750 3327
## delta_p:n5000 1.00 6067 3670
## causal_support:vistext_FALSE 1.00 1634 2196
## delta_p:vistext_FALSE 1.00 6644 3673
## n500:vistext_FALSE 1.00 3878 3940
## n1000:vistext_FALSE 1.00 3937 3412
## n5000:vistext_FALSE 1.00 4246 3907
## causal_support:delta_p:n500 1.00 6397 3544
## causal_support:delta_p:n1000 1.00 5542 3711
## causal_support:delta_p:n5000 1.00 4917 4021
## causal_support:delta_p:vistext_FALSE 1.00 6762 3793
## causal_support:n500:vistext_FALSE 1.00 2054 2876
## causal_support:n1000:vistext_FALSE 1.00 1841 2638
## causal_support:n5000:vistext_FALSE 1.00 1657 2350
## delta_p:n500:vistext_FALSE 1.00 5692 3882
## delta_p:n1000:vistext_FALSE 1.00 8084 3613
## delta_p:n5000:vistext_FALSE 1.00 7685 3420
## causal_support:delta_p:n500:vistext_FALSE 1.00 6769 3631
## causal_support:delta_p:n1000:vistext_FALSE 1.00 6332 3928
## causal_support:delta_p:n5000:vistext_FALSE 1.00 5713 3886
## sigma_abscausal_support 1.00 3062 3742
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m6)
Posterior predictive
model_df %>%
group_by(n, vis, workerId) %>%
data_grid(
causal_support = quantile(model_df$causal_support, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 20))),
delta_p = quantile(model_df$delta_p, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 10)))) %>%
add_predicted_draws(m6, prediction = "lrr_rep", seed = 1234, n = 500) %>%
mutate(
# transform to probability units
response_A_rep = plogis(lrr_rep) * 100
) %>%
ggplot(aes(x = causal_support, y = response_A_rep, fill = vis)) +
stat_lineribbon(.width = c(0.95)) +
geom_point(data = model_df, aes(y = response_A), alpha = 0.5) +
scale_fill_brewer(type = "qual", palette = 2) +
scale_color_brewer(type = "qual", palette = 2) +
labs(subtitle = "Posterior predictive distribution") +
theme(panel.grid = element_blank()) +
facet_wrap(workerId ~ .)
Model comparison
loo(m1, m2, m3, m4, m5, m6)## Warning: Found 11 observations with a pareto_k > 0.7 in model 'm5'. With this
## many problematic observations, it may be more appropriate to use 'kfold' with
## argument 'K = 10' to perform 10-fold cross-validation rather than LOO.## Warning: Found 15 observations with a pareto_k > 0.7 in model 'm6'. With this
## many problematic observations, it may be more appropriate to use 'kfold' with
## argument 'K = 10' to perform 10-fold cross-validation rather than LOO.## Output of model 'm1':
##
## Computed from 5000 by 306 log-likelihood matrix
##
## Estimate SE
## elpd_loo -581.9 13.6
## p_loo 3.2 0.5
## looic 1163.7 27.2
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm2':
##
## Computed from 5000 by 306 log-likelihood matrix
##
## Estimate SE
## elpd_loo -575.3 14.1
## p_loo 12.1 2.0
## looic 1150.5 28.3
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 304 99.3% 250
## (0.5, 0.7] (ok) 2 0.7% 491
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm3':
##
## Computed from 5000 by 306 log-likelihood matrix
##
## Estimate SE
## elpd_loo -576.3 13.5
## p_loo 19.1 2.6
## looic 1152.6 27.0
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 303 99.0% 602
## (0.5, 0.7] (ok) 3 1.0% 214
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm4':
##
## Computed from 5000 by 306 log-likelihood matrix
##
## Estimate SE
## elpd_loo -576.9 13.5
## p_loo 21.2 3.1
## looic 1153.8 27.1
## ------
## Monte Carlo SE of elpd_loo is 0.1.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 304 99.3% 490
## (0.5, 0.7] (ok) 2 0.7% 196
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm5':
##
## Computed from 5000 by 306 log-likelihood matrix
##
## Estimate SE
## elpd_loo -511.3 16.1
## p_loo 57.5 5.5
## looic 1022.7 32.2
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 270 88.2% 1042
## (0.5, 0.7] (ok) 25 8.2% 125
## (0.7, 1] (bad) 11 3.6% 65
## (1, Inf) (very bad) 0 0.0% <NA>
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'm6':
##
## Computed from 5000 by 306 log-likelihood matrix
##
## Estimate SE
## elpd_loo -514.5 15.8
## p_loo 63.9 5.7
## looic 1029.0 31.6
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 258 84.3% 501
## (0.5, 0.7] (ok) 33 10.8% 207
## (0.7, 1] (bad) 14 4.6% 56
## (1, Inf) (very bad) 1 0.3% 49
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## m5 0.0 0.0
## m6 -3.1 2.7
## m2 -63.9 9.8
## m3 -64.9 9.8
## m4 -65.6 10.0
## m1 -70.5 11.1Effects of interest as linear in log odds slopes
Derive linear in log odds slopes
slopes_df <- model_df %>%
group_by(n, vis, workerId) %>%
data_grid(
causal_support = c(0, 1),
delta_p = quantile(model_df$delta_p, probs = plogis(seq(from = qlogis(0.001), to = qlogis(0.999), length.out = 20)))) %>%
add_fitted_draws(m5, value = "lrr_rep", seed = 1234, n = 500) %>%
compare_levels(lrr_rep, by = causal_support) %>%
rename(slope = lrr_rep)Effect of sample size
slopes_df %>%
group_by(n, .draw) %>% # group by predictors to keep
summarise(slope = weighted.mean(slope)) %>% # marginalize out delta_p and vis effects
ggplot(aes(x = slope, y = "")) +
stat_halfeyeh() +
theme_bw() +
facet_grid(n ~ .)## `summarise()` regrouping output by 'n' (override with `.groups` argument)
Interaction of visualization and sample size
slopes_df %>%
group_by(n, vis, .draw) %>% # group by predictors to keep
summarise(slope = weighted.mean(slope)) %>% # marginalize out delta_p effects
ggplot(aes(x = slope, y = vis)) +
stat_halfeyeh() +
theme_bw() +
facet_grid(n ~ .)## `summarise()` regrouping output by 'n', 'vis' (override with `.groups` argument)
Pattern of slopes across levels of delta_p
slopes_df %>%
group_by(delta_p, .draw) %>% # group by predictors to keep
summarise(slope = weighted.mean(slope)) %>% # marginalize out n and vis effects
ggplot(aes(x = slope, y = delta_p, group = .draw)) +
geom_line(alpha = 0.1) +
theme_bw()## `summarise()` regrouping output by 'delta_p' (override with `.groups` argument)
Interaction of visualization and delta_p
slopes_df %>%
group_by(delta_p, vis, .draw) %>% # group by predictors to keep
summarise(slope = weighted.mean(slope)) %>% # marginalize out n and vis effects
ggplot(aes(x = slope, y = delta_p, group = .draw)) +
geom_line(alpha = 0.1) +
theme_bw() +
facet_grid(. ~ vis)## `summarise()` regrouping output by 'delta_p', 'vis' (override with `.groups` argument)
Effect of visualization
slopes_df %>%
group_by(vis, .draw) %>% # group by predictors to keep
summarise(slope = weighted.mean(slope)) %>% # marginalize out delta_p effects
ggplot(aes(x = slope, y = vis)) +
stat_halfeyeh() +
theme_bw()## `summarise()` regrouping output by 'vis' (override with `.groups` argument)